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Abstract 

Background: Many problems in biomedicine and other areas of the life sciences can be characterized as 
control problems, with the goal of finding strategies to change a disease or otherwise undesirable state of a 
biological system into another, more desirable, state through an intervention, such as a drug or other 
therapeutic treatment. The identification of such strategies is typically based on a mathematical model of the 
process to be altered through targeted control inputs. This paper focuses on processes at the molecular level 
that determine the state of an individual cell, involving signaling or gene regulation. The mathematical model 
type considered is that of Boolean networks. The potential control targets can be represented by a set of nodes 
and edges that can be manipulated to produce a desired effect on the system. 

Results: This paper presents a method for the identification of potential intervention targets in Boolean 
molecular network models using algebraic techniques. The approach exploits an algebraic representation of 
Boolean networks to encode the control candidates in the network wiring diagram as the solutions of a system 
of polynomials equations, and then uses computational algebra techniques to find such controllers. The control 
methods in this paper are validated through the identification of combinatorial interventions in the signaling 
pathways of previously reported control targets in two well studied systems, a p53-mdm2 network and a blood 
T cell lymphocyte granular leukemia survival signaling network. Supplementary data is available online and our 
code in Macaulay2 and Matlab are available via http://www.ms.uky.edu/~dmu228/ControlAlg. 

Conclusions: This paper presents a novel method for the identification of intervention targets in Boolean 
network models. The results in this paper show that the proposed methods are useful and efficient for 
moderately large networks. 
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1 Background 

The dynamics of gene regulatory networks (GRNs) 
have been studied using different modeling frame¬ 
works, with the goal of building computational mod¬ 
els of GRNs to get new insights into important cellu¬ 
lar processes, e.g., the cell cycle, [1,2], oscillations in 
the p53-mdm2 system, [3-5], the phage-lambda sys¬ 
tem, [6-8], or the T cell large granular lymphocyte 
(T-LGL) leukemia network, [9,10]. A generally diffi¬ 
cult problem is to design control policies to achieve 
desired dynamics in GRNs. This is particularly impor¬ 
tant in the control of cancer cells, [5,11-14] and cell 
fate reprogramming, [15,16]. Thus, developing tools to 
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control mathematical models of GRNs are key to the 
design of experimental control policies. 

There is a rich theory for the control of contin¬ 
uous models such as systems of differential equa¬ 
tions, [17-20]. Discrete models such as Boolean net¬ 
works (BN) have been proposed to study GRNs. In a 
BN, the genes of a GRN are represented by a set of 
nodes that can take on only two possible states (ON 
or OFF). Time is discrete, and the state of a gene at 
the next time step is determined by a Boolean function 
over a subset of nodes of the BN. The dependence of a 
gene on the state of another gene is graphically repre¬ 
sented by a directed edge. BN models are widely used 
in molecular and systems biology to capture coarse¬ 
grained dynamics of a variety of regulatory networks 
and have been shown to provide a good approximation 
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of the dynamics of continuous processes [8,10,21-28]. 
However, control methods for discrete models are still 
in their infancy, compared to the theory for continuous 
models. 

In this work, we propose a framework to study the 
control of BNs. Therapeutic interventions are modeled 
by manipulating the wiring diagram of a BN accord¬ 
ingly. For example, a gene knockout is modeled by fix¬ 
ing the value of one node of the BN to OFF. Similarly, 
deleting directed edges of a BN models the action of 
a drug that inactivates the interaction between two 
gene products. The control problem consists of finding 
the sequence of actions (node and edge deletions) to 
get the BN out of an undesirable state, and drive it 
towards a desirable state. Undesirable states may rep¬ 
resent a disease condition of a cell such as, for exam¬ 
ple, a highly proliferative state of a cancer cell, and a 
desirable state may represent cell death. Thus, a ther¬ 
apeutic intervention could aim at inducing the death 
of aberrant tumor cells. 

Several approaches to address this problem have 
been used in recent years. For example, the optimal 
control techniques developed in [29-32] assume a set of 
control nodes to derive a control policy that minimizes 
the likelihood of transitioning into an undesirable state 
in a stochastic context. Other approaches for the iden¬ 
tification of intervention targets include the use of sta¬ 
ble motifs in the network to induce the system into a 
desirable state [33]. In [34], integer programming was 
used to find a set of nodes to control the states of 
BNs representing tumor and normal cells. Optimiza¬ 
tion techniques were used in [35] to determine possible 
node manipulations for BNs. There are also search al¬ 
gorithms based on genetic and greedy algorithms de¬ 
scribed in [36-38]. For continuous models, a related 
control approach is given in [19]. 

The idea behind our approach is to encode the prob¬ 
lem of finding control candidates as a problem of solv¬ 
ing a system of polynomial equations. Then we can use 
computational algebra techniques to solve the system. 
This approach takes advantage of the rich algorithmic 
theory of computer algebra (e.g. Grobner basis tech¬ 
niques) and the theoretical foundations of algebraic ge¬ 
ometry with software available for computations (e.g., 
Macaulay2 [39]). The output of our method is a set 
of candidate actions with the capacity to induce the 
GRN towards desirable states. The method also has 
the ability of identifying combinatorial interventions 
in the network which are sometimes more effective as 
will be shown in the results section. The algebraic view 
of discrete models has been previously used for the de¬ 
velopment of computational tools to determine the at¬ 
tractors of BNs [40-43] , and also for inferring network 
structure from dynamics [44-47]. 


2 Methods 

2.1 Boolean Networks 

A Boolean network is a dynamical system that is dis¬ 
crete in time as well as in variable states. More for¬ 
mally, consider a collection xi,..., a;„ of variables that 
take values in the binary set {0,1}. Then a Boolean 
network in the variables xi,..., x„ is a function 

F = (/!,...,/„): {0,1}" ^{0,ir 

where each coordinate function fi : {0,1}" —>■ {0,1} is 
a Boolean function on a subset of {xi,..., x„} which 
represents how the future value of the *-th variable 
depends on the present values of the other variables. 

Example 2.1 For concreteness, we illustrate the defi¬ 
nitions using the following toy network 

/i = ^X3A^X5, /2 = ^a;iVX4, f3=-'X2, 

fi = X3, /5 = -1X4, 

where A, V, and ^ are the AND, OR, and NOT log¬ 
ical operators, respectively. In the context of model¬ 
ing biological systems, A corresponds to activation by 
the combination of regulators (all regulators are nec¬ 
essary), V corresponds to independent activation (one 
regulator is sufficient), and ^ corresponds to negative 
regulation. 

Given a Boolean network F = (/i,...,/„), a di¬ 
rected graph yy with n nodes xi,..., x„ is associated 
with F. There is a directed edge in W from Xj to Xi 
if Xj appears in fi. Notice that the presence of the 
interaction Xj —>■ Xi implies that the regulatory func¬ 
tion fi depends on Xj , say fi (x^j,..., Xj,..., x*,^) with 
Xj € {xfcj,..., Xk^}. In the context of a molecular net¬ 
work model, yy represents the wiring diagram of the 
network. 

Example 2.1 (cont.) The wiring diagram of the toy 
network is shown in Figure lA. 

If the set {0,1} is given the structure of a finite 
field with standard addition and multiplication, F 2 = 
{0,1}, then the functions fi : ¥2 ^ ¥2 can be rep¬ 
resented as polynomials over F 2 , and the dynamical 
system F = (/i,..., /„) : FJ —>■ F 2 becomes a polyno¬ 
mial dynamical system, see [41], which gives access to 
a range of mathematical tools for the analysis of F. 

Example 2.1 (cont.) The rules to transform Boolean 
functions into polynomials in F 2 [xi,..., x„] are given 
hy a A b = ab, a V b — a + b + ab, and = 1 -I- a, 
where the operations are computed modulo 2. For the 
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A 

©■— 

molecular network 
(u=control parameters) 
X(t+l)=jr(x(t),u) 


B 



Desired attractor 
landscape: L* 


desired 
attractor 

Attractor landscape Lp 
(no control, F(x)=,7^(x,0)) 

I Polynomial representation 
of molecular network 
Ji^(x,u)GF2[x,u] 

Identification of control 
targets: solve Ljp-(x,u)=L* 

algebra 

=A, B, C, D, 


system will get stabilized are of particular importance. 
These special points of the state space are called at¬ 
tractors of a Boolean network and these may include 
steady states (fixed points), where F(a;) = x, and cy¬ 
cles, where F^'(x) = x for some integer A: > 1. At¬ 
tractors in Boolean network modeling might represent 
cell types [48] or cellular states such as apoptosis, pro¬ 
liferation, or cell senescence [49,50]. Identifying the 
attractors of a system is an important step towards 
the control of that system and can be done using tools 
from computational algebra [40,41]. 

Example 2.1 (cont.) The steady states of the toy 
network are found by solving the system of equations 
/i = Xi, f = 1,..., 5. This means we want to find the 






roots of gi = 0, where gt = fi — Xi. This gives 

the 

% 

■ 

system of equations 


u=A 


u = B 


<?i = 1 + 2:3 -l- X 5 -|- X 3 X 5 -|- Xi = 0, 





‘AM’ 

gl 2 = 1 -1- Xi -1- X 1 X 4 -1- X 2 = 0, 

^3 = 1 -1- X 2 -1- X 3 = 0, 
g4 = X3 + XA = 0, 

(1) 

u=C 


u = D 


35 = 1 -1- X 4 -1- X 5 = 0. 



Figure 1 Description of the algebraic approach of 
identification of control targets. A. A BN model of a molecular 
network. The control variables are the entries of u. B. In the 
absence of control policies {u = 0), the attractor landscape 
{Lp) can have undesirable attractors. C. The goal is to 
choose control values that give a desired attractor landscape, 
L*. D. To use the algebraic approach we first find the 
polynomial representation of the BN (see Section 2.2). E. The 
next step is to set up the desired attractor landscape as a 
system of equations that the BN has to satisfy, = L* 

(see Section 2.3). F. Solving the equation 

will provide the control values to achieve the desired landscape 
(see Section 2.4). This approach not only finds individual 
control policies {u = A, single node; u = B, single edge), but 
also combinatorial control policies {u = C, two nodes; u = D, 
two edges). In a combinatorial control policy, the desired 
attractor landscape is achieved by the combination of two or 
more entries of u. 


Note that, since the system is not linear, we cannot 
make use of tools such as Gaussian elimination. Com¬ 
putational algebra allows us to solve such systems by 
encoding the solutions as an algebraic object called an 
ideal of polynomials, I = { 91 , 92 ^ 93 ^ 94 -, 95)1 and then 
finding an equivalent, but simpler representation. More 
precisely, for the system in Equation 1 we can find its 
Grdbner basis [51] and obtain 

I = (91,92,93,94,95) = {X 5 + ^,X 4 ,X 3 ,X 2 + l,Xi). 

This means that the original system has the same 
solutions as the simpler system 

xi = 0, X 2 -I- 1 = 0, X 3 = 0, X 4 = 0, X 5 -f 1 = 0, 


toy network we obtain 

/l = (1 -I- X3)(l -I- X5) = 1 -I- X3 -I- X5 -I- X3X5, 

/2 = (1 + Xi) -I- X4 -I- (1 -I- Xi)xA = 1 -I- Xi -I- X1X4, 
/3 = l+X2, 

/4 = 2 : 3 , 
f3 = l+ X 4 . 

The dynamical properties of a Boolean network are 
given by the difference equation x{t -I- 1) = F(x(t)); 
that is, the dynamics is generated by iteration of F. 
More precisely, the dynamics of F is given by the state 
space graph S, defined as the graph with vertices in 
F 2 = { 0 , 1 }", which has an edge from x € ¥2 to y € ¥2 
if and only if y = F(x). The states x € F 2 where the 


which is easily solved to obtain x = 01001 as a 
steady state of the toy network (parentheses omitted 
for brevity). 

2.2 Control Actions: edge and node manipulations 
This paper considers two types of control actions: 1. 
Deletion of edges and 2. Deletion (or constant expres¬ 
sion) of nodes. The motivation for considering these 
actions is that they represent the common interven¬ 
tions that prevent a regulation from happening. For 
instance, an edge deletion can be achieved by the use 
of therapeutic drugs that target specific gene inter¬ 
actions and node deletions represent the blocking of 
effects of products of genes associated to these nodes; 
see [5]. 








Murrugarra et al. 


Page 4 of 12 


A schematic of our approach is given in Figure 1 
where the dynamics of a Boolean network can be ma¬ 
nipulated by a set of controls consisting of deletion 
or constant expressions of edges and nodes. Below we 
define and explain all the steps in more detail. 

A Boolean network with control is given by a func¬ 
tion : F 2 X 17 —f F 2 , where 17 is a set that denotes all 
possible control inputs (Figure lA). Given a control u 
in 17, the dynamics are given by x(l -|- 1) = T{x{t),u). 

We consider a BN F = (/i,..., /„) : F 2 —f F 2 and 
show how to encode edge and node controls by : 
F 2 X [/ —f F 2 , such that T{x,0) = F(a;). That is, the 
case of no control coincides with the original BN. We 
remark that these control types can be combined, but 
for clarity we present them separately. 

Definition 2.2 (Edge Control) Consider the edge 
Xi —f Xj in the wiring diagram W. The function 

Tj{x,u^j) := fj{xi,{uij + 1 ) 2 :*, ■..,Xn) ( 2 ) 

encodes the control of the edge Xi —f Xj , since for each 
possible value of uij S F 2 we have the following control 
settings: 

• If Uij = 0, Tj{x, 0) = /j(xi,... ,Xi,... ,Xn)- That 
is, the control is not active. 

• If Ui^j I, J'j (x,l) ^j(xi,...,Xj 0,..., Xji') . 

In this case, the control is active, and the action 
represents the removal of the edge Xi —>■ Xj. 

Definition 2.2 can easily be extended for the control of 
many edges, so that we obtain 7^ : F 2 xF^ —f F 2 , where 
e is the number of edges in the wiring diagram. Each 
coordinate, Uij, of u in T{x,u) encodes the control of 
an edge Xi —f Xj. 

Example 2.1 (cont.) Incorporating edge control in 
the toy network results in 

J -1 = 1 + (lt3_i -P l)x3 -p (rt 5 p -P l)X5-p 
(U 3 ,l -P 1 )X 3 (U 5,1 -P 1 )X 5 , 

7^2 = 1 -P (m 1,2 + l)xi -p (mi ,2 + l)xi(w4,2 + 1)X4, 
^3 = 1 -p (rt 2,3 + 1 )X 2 , 

-^4 = (W 3.4 + 1)X3, 

T's = 1 -p (U 4,5 -p 1 )X 4 . 

Definition 2.3 (Node Control) Consider the node 
Xi in the wiring diagram W. The function 

Fj (x, u ~, u+) := {u~ +ut + l)fj (x) -P ut (3) 

encodes the control (knock-out or constant expres¬ 
sion) of the node Xi, since for each possible value of 
{u~ ,uf) S F 2 we have the following control settings: 

• For u~ = 0,ull = 0, 7(,(x,0,0) = fj{x). That is, 
the control is not active. 


• For u~ = l,ul = 0, Fj{x, 1,0) = 0. This action 
represents the knock out of the node Xj. 

• For u~ = 0,uf = 1, 7^j(x, 0,1) = 1. This action 
represents the constant expression of the node Xj. 

• Forrt” = l,u+ = 1, J)(x, 1,1) =/^(xti,... ,Xt„)-P 
1. This action changes the Boolean function to its 
negative value and might not be a relevant case 
of control. 

Example 2.1 (cont.) Incorporating node control in 
the toy network results in 

Fi = {u^ +Ui + 1)(1 -P X 3 -P X 5 -P X 3 X 5 ) -P uf, 

= {U 2 +U 2 + 1)(1 -P Xi -P X 1 X 4 ) -P U 2 , 

^3 = (^3 + ■^^ + + * 2 ) + ) 

J 4 = (^4 -P -P 1)X3 -P U 4 , 

^5 = {^5 + -P 1)(1 -P X 4 ) -P 

2.3 Control targets in Boolean networks. 

We consider a BN with control F : ¥2 x U —>■ ¥2 , and 
denote by F the BN with no control (7^(x, 0) = F(x)). 
We remark that in each case of interest, both edge and 
node control could be analyzed simultaneously. 

2.3.1 Generating new steady states 
Suppose that yg = {yoi, ■ ■ ■ ,yon) S F 2 is a desirable 
cell state (for instance, it could represent the state of 
cell senescence) but is not a fixed point, i.e., F(yg) ^ 
yg. The problem is then to choose a control u such that 
■^(yo, u) = yo- We now show how this can be achieved 
in our framework. 

After encoding our BN with control as a polynomial 
system Fj{x,u) G F 2 [x, 4 t] (see Section 2.2), we con¬ 
sider the system of polynomial equations in the u pa¬ 
rameters: 

(Yo: u) - Voj = 0, j = 1,..., m. (4) 

Example 2.1 (cont.) Here we consider a toy network 
and assume we are interested in controlling edges to 
make yg = 01111 a steady state. For simplicity we 
only consider control using the edges Xi —)■ X 2 ,X 4 —)■ 
X 2 , X 2 —>■ X 3 , X 4 —)■ X 5 . The network is 

7^1 = 1 -P X 3 -P X 5 -P X 3 X 5 , 

F2 = 1 -p {ui ^2 + 1)3^1 + (^ 1,2 + l)xi(u 4_2 + 1 )X 4 , 
F 3 = 1 -p (1*2,3 + 1)3:2: 

Fa = X 3 , 

J's = 1 -P (U4,5 -P 1)X4. 

Evaluating at yg =01111 we obtain 

Fl =0, F 2 = 1, F 3 = 4*2,3: -^4 = 1: F 5 = 4 * 4 , 5 . 
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Then, 01111 will be a steady state if and only if Ti = 
Yqj for i = 1,..., 5. This gives the following solution 
U 2,3 = 1, U 4,5 = 1. Thus, in this case there is a unique 
control, U2,3 = ^4,5 = 1, that guarantees yg = 01111 
is a steady state (this control policy is illustrated in 
Figure IF, u = D). In general, the equations can be 
more complex and nonlinear, so computational algebra 
is needed. See the results section for applications to 
more complex models. 

Example 2.1 (cont.) We now show how the problem 
of creating new steady states can also be solved us¬ 
ing node controls. We again use the toy network, and 
assume that we want to make yg = 11110 a steady 
state. For simplicity we only consider control of nodes 
xi (constant expression), X3 (knock-out and constant 
expression), and node X4 (constant expression). The 
network is 

-^1 = (■*^1' + 1)(1 + 2^3 -f Xg -|- X3X3) + uf , 

7^2 = 1 + 2;i -I- X 1 X 4 , 

^3 = ('*23 -|- -|- 1)(1 -|- X2) + , 

J4 = (uj; + 1 )X 3 + uj, 

E 3 = 1 -|- X4. 

Evaluating at yg = 11110 we obtain 

J'l=U^, 7^2 = 1, 7^3 = Mg, 7^4 = 1, 7^5 = 0. 

Then, 11110 will be a steady state if and only if = 
1 and = 1. That is, neither control by itself is 
sufficient, but together they create the steady state 
(this control policy is illustrated in Figure IF, u = C). 

3.3.3 Destroying existing steady states, or, in 
general, blocking transitions 
Suppose that Xg S F 2 is an undesirable steady state of 
F(a:), that is, F(xg) = Xg (for instance, it could rep¬ 
resent a tumor proliferative cell state that needs to be 
avoided). The goal here is to find a set of controls such 
that 7^(xg,u) ^ Xg. More generally, we may want to 
avoid a transition between two states Xg and zg. That 
is, we want to find controls such that 7^(xg,u) ^ zg. 
To solve this problem consider the following equation, 

(7'l(Xo,u)-Zoi-bl) . . . (J'„(xg,u)-Zg„-bl) = 0. (5) 

Equation 5 defines a polynomial equation in the u pa¬ 
rameters. It can be shown that 7^(xg,u) ^ Zg if and 
only if Equation 5 is satisfied. 

Example 2.1 (cont.) Here we focus on finding edges 
to block the transition from xg = 00101 to F(xg) = 
01111. For simplicity we only consider control using 


the edges X3 —)■ a;i,a;5 —>■ X4,X2 —>■ X3,X3 —)■ a;4. The 
network is 

7^1 = 1 -I- (u3,i -I- l)x3 -I- (wsq -I- l)a;5-|- 

(m3.i + l)a;3(M5.i + 1 ) 2 : 5 , 

^3 = 1 + {U 2,3 + 1)2:2, 

7^4 = (m3,4 + 1 ) 2 : 3 , 

^3 = 1 + X4. 

Evaluating at Xg = 00101 we obtain 

•^1 = M3pU5q,7^2 = 1,7^3 = 1,7^4 = M34,7^5 = 1. 

Then, Eq. 5 becomes (m3,iM5,i- b 1 )(123,4-|- 1 ) = 0 which 
has two solutions: 223,4 = 1 and 223,1 = 225,1 = 1 (first 
control policy illustrated in Figure IF, u = B). The 
second solution is what we refer to as a combinatorial 
control policy. Neither 223,1 = 1 or 225,1 = 1 is sufficient, 
but combined they block the transition. In general, 
the equations can be more complex and nonlinear, so 
computational algebra is needed. 

3.3.3 Blocking regions in the state space 
We now consider the case where we want the dynamics 
to avoid certain regions. For example, if a particular 
value of a variable, Xk = a G V2 triggers an undesirable 
pathway, or is the signature of an abnormal cell, then 
we want all steady states of the system to satisfy Xk ^ 
a. In this case, we consider the systems of equations 

Ej{x,u) - Xj = 0 ,j = l,...,m, 

n 

Xk — a = (). 

Note that, in contrast to Sections 2 . 3.1 and 2 . 3 . 2 , we 
are now using variables for x instead of specific values. 
Since the steady states with Xk = a are to be avoided, 
we want to find controls 22 for which Equation 6 has 
no solution. 

Example 2.1 (cont.) Here we consider the toy net¬ 
work and focus on controlling nodes to avoid regions 
of the form X3 = 0 . Eor simplicity we only consider 
control using the nodes X2 (knock-out), X3 (constant 
expression), X4 (constant expression). The network is 

Fi = 1+X3+X3+ 2:3X5, 

7^2 = {U 2 + 1)(1 +2:1-1- X1X4), 

+3 = {ut + 1)(1 + 2 : 2 ) + M 3 , 

E4 = (M4 -b 1)X3 -b m|, 

+■5 = 1 -b X4. 
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Then, Eq. 6 becomes 

1 + ^3 + Xs + 0:3X5 + Xi = 0, 

{U2 + 1)(1 + Xi + X1X4) + X2 = 0 , 

+ 1)(1 + X2) + + X3 = 0, 

(U4 + 1)X3 + U4 + X4 = 0 , 

1 + X4 + X5 = 0, 

X3 = 0 . 

In contrast with the previous examples, this system 
of equations cannot be analyzed by hand. In Sec¬ 
tion 2.4 we will show how computational algebra gives 
an equivalent, but simpler, system of equations. 

2.4 Identifying control targets 

In each case of Section 2.3 we obtained a system of 
equations (or a single equation) that we need to solve 
to find the appropriate controls. This can be done us¬ 
ing computational algebra tools. For instance, we can 
compute the Grobner basis of the ideal associated to 
Equation 4, see [51], 

1= {■^i{yo,u) - yoi, ■ ■ ■ ,J^n{yo,u) - yon) ■ ( 7 ) 


Results 

We test our control methods using two published mod¬ 
els where potential intervention targets were identified 
along with experimental validations. In the first ex¬ 
ample we focus on control with edge deletions while 
for the second we use control with node deletions and 
constant expressions. 

Example 2.4 P53-mdm2 network. A Boolean net¬ 
work model for the signaling response of DNA damage 
in a p53 network was developed in [5]; the wiring di¬ 
agram of this model is reproduced in Figure 2. The 
tumor suppressor protein p53 can trigger either cell 
cycle arrest or apoptosis in response to DNA dam¬ 
age in normal cells. The authors of [5] extended their 
analysis to a breast cancer cell line, MCF7, with the 
purpose of identifying potential therapeutic interven¬ 
tions in the network for the cancer cell. Thus, in this 
example we will focus on the cancer cell model where 
PTEN and pl4ARf are always inactive (fixed to zero) 
and cyclinG is always active (fixed to 1), see Table 5 
of [5] . This system can be represented as a discrete dy¬ 
namical system F = (fi,, /ig) : F 2 ® —?► F 2 ® with 16 
nodes and 50 edges. We represent the nodes by 


Example 2.1 (cont.) Now we continue the previous 
example where the goal was to avoid regions of the 
form X3 = 0 . We encode the system of equations as 
I = (1 -I-X3 -I-X5 -I-X3X5 -bXi, {U2 + 1)(1 -tXi -I-X1X4) -I- 
X2, + 1)(1 -I- X2) + + X3, (^4 -I- 1)X3 -I- M4 -I- X4, 

1 -I- X4 -|- X5, X3). Using computational algebra we find 
a Grobner basis of this ideal: 

I = {uf, U2 ,X5 -I-U4 -b 1 ,X 4 -!-«;{■, X3,X2 -b l,Xi -b M4 ) . 

Thus, the original system of equations has the same 
solutions as the system 

= 0 , U2 = 0, X5 -b -b 1 = 0 , X4 -b MJ = 0 , 

2:3 = 0, X 2 -b 1 = 0 , xi -b W4 = 0. 

In order to avoid regions of the form X3 = 0 , we 
need to find parameters for which the system has no 
solutions. This is guaranteed if any of the polynomials 
is equal to 1 . Namely, we select the equations that only 
have the control parameters = 0 , = 0 . If we use 

u]]' = 1 or it2 = 1, the system is guaranteed to have no 
solutions. Thus, = 1 or = 1 independently are 
sufficient to achieve the control of the system (second 
control policy illustrated in Figure IF, u = A). 

As we will show in Examples 2 . 4 - 2 . 5 , the computa¬ 
tion of a Grobner basis allows us to read out all con¬ 
trols for which the system of equations has a solution. 
Furthermore, our algebraic approach can detect com¬ 
binatorial control actions (control by the synergistic 
combination of more than one action). 


xi = ATM, 

X 3 = Mdm2, 

X 5 = Wipl, 

X 7 = PTEN, 
Xg = AKT, 

xii = Rh, 

a:i3 = pUARJ, 
X 15 = Bax, 


X 2 = p53, 

X 4 = MdmX, 
Xg = cyclinG, 
xs = p21, 
xio = cyclinE, 
X 12 = E2E1, 
Xi 4 = Bcl2, 
Xig = caspase. 


( 8 ) 


The polynomial functions for this network are listed 
in the supplementary material. We remark that the 
original model in [5] considers threshold functions for 
all the regulatory rules and this type of functions might 
introduce self-loops in the wiring diagram. That is, 
when we translated the threshold rule for Xi into a 
polynomial fi, the function fi might depend on x^. 
Notice that the self-loops were omitted in Figure 2 to 
be consistent with the original model. 

For this model, in the presence of DNA damage, the 
system has a single limit cycle representing the state 
of cell cycle arrest, where p53 and p21 are oscillating; 
see Figure 3. 

Suppose that we want to find out which set of edges 
one can manipulate in order to destroy this limit cycle 
and to lead the system into a different attractor, one 
that represents a desirable cell state. Let us start by 
considering a desirable state that represents cell death, 


yo = (l,l,l,0,0,l,0,l,l,0,0,l,0,0,l,l). 


(9) 


Murrugarra et al. 


Page 7 of 12 



Figure 2 The p53-mdm2 network adapted from [5]. Arrows in 
green represent activation while hammerhead arrows (in red) 
represent inhibition. Self loops were omitted, see text in 
Example 2.4 for an explanation. For the cancer cell model, 
PTEN and pl4ARf are always inactive (fixed to zero) and 
cyclinG is always active (fixed to 1). 


x7=0011010010010100 

xl=1010010011010100 x6=0011110110010100 

I t 

x2=1000010011010100 x5=0111110110110100 

t f 

x3 = 1101010011110100^>x4=1101110111110100 

Figure 3 States of limit cycle representing cell cycle arrest in 
the p53 model. The order of the vector entries follows the 
indexing in Equation 8. The dashed edge represents the 
transition target to destroy the limit cycle. 


where Xig = caspase is active. In order to make yg a 
steady state with the approach described in the Meth¬ 
ods section, we form the following system of polyno¬ 
mial equations, 

-^*(yo>Ai) = = 1, ■ • ■, 16. (10) 

The solutions for the system of equations 10 are given 
by the nonzero generators of the ideal associated to 
the system 10, given in Equation 11, 

{^^2,5 + Ij + 1 ), Us,8(^3,8 + 1 ); U2,2{U3,2 + 1 ), 

Ul,2iu^,2 + 1)) Mi2,16'^16,16('^8,16 + 1)} 

( 11 ) 

There are 945 solutions for Equation 11. Each solution 
gives a controller that guarantees the desired steady 
state, but each controller may have a different impact 


Table 1 Difference of impact in the combinatorial action of edge 
deletions. Control edges that increase the basin of attraction of 
cell death represented by yg in Equation 9. There are 
2^® = 65536 possible states. The number in parentheses is the 
ratio between the basin size and the total number of states. 


Controllers applied 

Ref. 

Basin size of Yo 

mdm2 p53 
p53 Wipl 

[5] 

35581 (54.29%) 

p53 —>■ Wipl 

Mdm2 —> p21 
Mdm2 —>■ p53 
p21 —)■ Caspase 

A control set that 
forces Yg to be 
a fixed point, from 
Equation 11. 

39856 (60.82%) 

mdm2 —>■ p53 
p53 —> Wipl 
mdm2 — p21 
p21 —> Caspase 
ATM Rb 
mdm2^ Rb 
mdmx p53 

Rb E2F1 

Bcl2 Bax 

A control set to make 
Yg a fixed point 
and for blocking 
the transition in red 
at Figure 3. 

65536 (100%) 


on the dynamics of the network. Table 1 (middle row) 
shows one of the controllers from Equation 11. In Ta¬ 
ble SI of the supplementary material, we list the top 
ten control combinations from Equation 11, that is, 
those that give the largest basin of yg as well as the ten 
sets that give the smallest basin for yg; see Table S2. 
For each solution in tables S1-S2, we crossed-out the 
edges that become nonessential after the other con¬ 
trollers in the set have been applied; see the discussion 
for more about nonessential edges. 

Furthermore, we can aim to destroy the limit cycle 
in Figure 3. We target one of the transitions within 
this limit cycle, let us take the transition Xg —f Xg, see 
dashed transition in Figure 3. Blocking this transition 
gives additional controllers that help to stabilize the 
system at the desired fixed point; see the last row of 
Table 1. 

The deletion of the edges mdm2 —f p53 and p53 —f 
Wipl were identified in [5] by deleting one edge at 
a time and checking if the modified system had the 
desired attractor. [5] reported that the combinatorial 
action of these controllers increased the basin of at¬ 
traction of the desired fixed point to more than 50%, 
which was validated experimentally. However, doing 
this type of search for finding all possible combinations 
of controls is infeasible. Since for each edge we have 2 
possible actions (control or no control), and there are 
50 edges, then there are 2®° networks to be analyzed in 
total. In contrast, the computational algebra approach 
of this paper allows to obtain all combinations of edges 
that guarantee the desired steady state of the system 
in one process. 
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Example 2.5 T-LGL network. A Boolean network 
model for the blood cancer large granular lymphocyte 
(T-LGL) leukemia was developed in [9] and further 
analyzed in [10,33]. T-LGL leukemia is characterized 
by escaping cell death through abnormal mechanisms, 
which are insensitive to Fas-induced apoptosis, [9]. 
This network has 60 nodes but was reduced to a sub¬ 
network of 16 nodes for steady state analysis, see Fig¬ 
ure 4. Here we use the 16-nodes network to identify 
potential control targets. 



Figure 4 Reduced T-LGL network adapted from [10]. Arrows 
in green represent activation while hammerhead arrows (in 
red) represent inhibition. All the negative edges from 
Apoptosis were omitted for simplicity. 


We represent the nodes by 
xi =CREB, X2 =IENG, 

X4 =GPCR, X5 =SMAD, 

X 7 =sFas, 
xio =Caspase, 

Xi3 =IAP, 

Xi6 = Apoptosis. 

The polynomial rules for this system are listed in 
Equation 12, 


X 3 =P 2 , 
xq =Fas, 
xs =Geramide, xg =DISC, 
xii =FLIP, xi2 =BID, 

Xi4 =MGL1, xi5 =S1P, 


fi — 2 ^ 2(1 + Tie), /2 — (xs -I- l)(x3 -|-1)(1 -I- xie), 

/a = (a;2a;3 + X 2 + X3)(l -I- xie), 

/4 = 2^15(1 + a;i 6 ), fb = 2:4(1 + 2:16), 

/e = ( 2:7 + 1)(1 + 2 ;i 6 ), /? = 2 : 15(1 + xie), 

/s = ( 2:15 + l)x 6 (l -bxie), 

/g = ( 2 : 6 X 3 X 11 -b XgXg -b XeXii -b Xe -b 2 : 8)(1 + Xie), 

/lO = (X 9 X 12 X 13 -b X 9 X 12 -b X 12 X 13 -b Xg -b Xi 2 )(l + Xie), 

/ll = ( 2:9 -b 1)(1 -b Xie), /12 = ( 2:14 + 1)(1 + 2 ;i 6 ), 

/l 3 = ( 2:12 + 1)(1 + 2 ;i 6 ), /i 4 = ( 2:9 -b 1)(1 + Xie), 

/l 5 = ( 2:3 + 1)(1 + 2 : 16 ), 

/16 = a:ioXi 6 + xio + xie- 


( 12 ) 

This system has three steady states, one that repre¬ 
sents the normal state, xg = 0000000000000001, where 
Apoptosis is ON, and 0001101000101110, 
0011101000101110, the disease states in which Gaspase 
and Apoptosis are OFF. 

To find the potential node deletions (or constant ex¬ 
pressions) that will block the disease states we use Eq. 
6 , 


-b u)*" -b l)/l (x) -b u)*" — Xi = 0 


(13) 

(■Uie + ufg -b l)/i6(x) -b - xie = 0 
xie = 0,xio = 0 


Equation 13 encodes all parameters for which there is 
a steady state where Gaspase and Apoptosis are OFF. 
Thus, we are interested in the parameters for which 
this system of equations has no solution. Since for each 
node we have 3 possible actions (no control, node dele¬ 
tion, constant expression), there are 3^® networks to be 
analyzed in total. Thus, even for this small network, 
exhaustive search is computationally challenging. 

On the other hand, computational algebra allows us 
to obtain the parameter combinations that guarantee 
that the disease states are not fixed points of the sys¬ 
tem (see the supplementary material for details). The 
parameter combinations (enclosed in brackets, where 
entries not shown are equal to zero) are 

= 1 }> {wg = 1 }. {wfo = {'^12 = 1 }. = 1 }> 

= 1}> {“6 = “n = 1}> {«7 = = !}• 

(14) 

Thus we obtain that the constant expression of Ce- 
ramide, DISC, Gaspase, or BID, or the deletion of 
MGLl or SIP, will guarantee that the disease states 
are not steady states of the system (Table 2). These 
controls could also be found by trying one control at 
a time as was done in [10]. Importantly, our compu¬ 
tational algebra approach shows that there are two 
additional control policies that consist of a combina¬ 
tion of different controls. It can be shown that neither 
Fas, sFas, FLIP individually can eliminate the disease 
states, but the deletion of FLIP combined with the 
constant expression of Fas or the deletion of sFas will 
work (Table 2). Furthermore, those are the only com¬ 
binations that guarantee the disease states will not be 
attractors. 

We note that, in the worst case, computing the 
Grobner basis for a system of polynomial equations 
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has doubly exponential complexity in the number of 
solutions. However, for the type of networks discussed 
in this paper, namely, biological networks where most 
of the nodes are regulated by only a small subset of 
the other nodes, Grobner bases can be computed in a 
reasonable time, see [42], 


Table 2 Control nodes for the reduced T-LGL network. The last 
two rows represent combinatorial actions of two nodes. All 
attractors are steady states, and the basin sizes include the steady 
states themselves. Notice that node xiq = Apoptosis is a 
conceptual node in this model, thus it is not a relevant solution 
for network control. 


Solution 

Control targets 

Attractor 

Basin size 


Ceramide—ON 

0000000100000001 

100% 

< = 1 

D/SC=ON 

0000000010000001 

100% 

■^10 = 1 

Caspase=ON 

0000000001000001 

100% 

'*^12 = 1 

e/D=ON 

0000000000010001 

100% 

Uii = 1 

MCLl=OFF 

0000000000000001 

100% 

^15 = 1 

SIP^OFF 

0000000000000001 

100% 

< = 1 
Mil = 1 

Fas^ON 

FLIP^OFF 

0000010000000001 

100% 

= 1 
Wn = 1 

sFas=OFF 

FLIP=OFF 

0000000000000001 

100% 


Discussion 

The design of control policies for gene regulatory net¬ 
works is an important challenge in systems biology. 
The method described here exploits the interplay be¬ 
tween the structure and the dynamics of the network to 
identify potential control interventions that will drive 
the system towards desired dynamics. The formula¬ 
tion of the problem as that of finding all solutions to 
a system of polynomial equations provides an alterna¬ 
tive to exhaustive search of all possible combinations 
of interventions, which often is not feasible. 

One shortcoming of the method is that, in its cur¬ 
rent form, it requires the Boolean network model to 
be updated synchronously, with deterministic dynam¬ 
ics. While steady states of the network do not de¬ 
pend on the update order used, general limit cycles 
and attractor basins do, however. Thus, some of the 
methods described here might not be applicable in a 
stochastic setting. For instance, if the dynamics is gen¬ 
erated using the asynchronous update or more gener¬ 
ally using the stochastic settings in [8,52-54], then 
encoding the controllers for blocking transitions will 
need to have a different setup than the one described 
here. However, the method for producing a new steady 
state is still valid for all variants of stochastic updates 
because the system will maintain the steady state. 
Another shortcoming is that the control methods de¬ 
scribed in this paper were developed for Boolean net¬ 
works only. However, many published discrete models 
of GRNs are multistate. These methods could poten¬ 
tially be extended to a more general setting, where the 
network variables might attain more than two states 
[41,54-57]. We remark that it is possible to map a mul¬ 
tistate model into a Boolean model (see [58]) where our 
methods can be applied and then it would be possible 
to recover a control set for the multistate model. 
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An important challenge in the process of identifi- 
cation of control targets in a network is to develop 
methods that can identify controllers that can guaran¬ 
tee global reachability of a desired steady state. This 
problem is not addressed in this paper, and remains 
to be studied in the future. For instance, the control 
strategies do not give any information about the basin 
size of a fixed point generated by the methods of this 
paper. However, we remark that some algebraic meth¬ 
ods allow to estimate the change in the basin size af¬ 
ter an edge deletion, see [59]. Nonetheless, the control 
targets identified by the algebraic techniques described 
here could be used for further analysis of the system, 
such as for studying reachability [60], or for designing 
optimal control policies in a stochastic setting [29-32] . 

Finally, it is worth pointing out that the methods 
of this paper might produce a large number of control 
strategies, which all give the desired result, but many 
of these might be biologically meaningless or infeasi¬ 
ble as actual interventions. Eliminating those might 
be challenging. For instance, some solution sets might 
contain nonessential [61] or nonfunctional [62] edges. 
That is, an edge could become nonessential after the 
other controllers in the set have been applied. In Ta¬ 
ble SI of the supplementary material, we list the top 
ten control combinations for Example 2.4 where we 
crossed-out the edges that become nonessential after 
the other controllers have been applied. Moreover, one 
can group the solution sets by considering only mini¬ 
mal sets where all the control edges are functional. For 
instance, we can group the first four solutions of Ta¬ 
ble SI of the supplementary material into one group 
with the minimal representative set given in the mid¬ 
dle row of Table 1 

Conclusions 

This paper presents a novel approach to the 
identification of potential interventions in Boolean 
molecular networks. The methods use the theoretical 
foundations of algebraic geometry to encode the 
structure of a network by a set of polynomials and 
then, with the use of computer algebra techniques, 
find a set of nodes and edges to perform interventions 
in silica. The methods were validated using two 
published models where dynamic network 
interventions were identified, the p53-mdm2 system 
and the T-LGL leukemia model. It was shown that 
the methods in this paper can identify the controllers 
that were already reported and also find new 
potential targets. Some of these new control targets 
are combinatorial in nature and might result in more 
efficient strategies as was shown in the results section 
using the change in the basin size of the system as a 
measure of efficiency. 
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